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We consider control of reaction and population systems by deterministically imposed transitions 
between the states with diflterent numbers of particles or individuals. Even where the imposed tran- 
sitions are significantly less frequent than spontaneous transitions, they can exponentially strongly 
modify the rates of rare events, including switching between metastable states or population extinc- 
tion. We also study optimal control of rare events, and specifically, optimal control of disease extinc- 
tion for a limited vaccine supply. A comparison is made with control of rare events by modulating 
the rates of elementary transitions rather than imposing transitions. It is found that, unexpectedly, 
for the same mean control parameters, controlling the transitions rates can be more efficient. 
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I. INTRODUCTION 

Many features of biophysical systems and population 
dynamics are related to the randomness of the underly- 
ing processes, such as elementary molecular reactions, or 
birth and death, or infection and recovery. The random- 
ness is important where the number of involved molecules 
or individuals, even though it is large, is not very large, 
as is often the case in gene expression, for example [1]. 
In such mesoscopic domain fluctuations are small on av- 
erage. Then, if the system is in stationary conditions, for 
much of the time its molecular composition or the popu- 
lation size experience small variations about the values in 
dynamically stable states, i.e., the stable states that the 
system vifould have with no fluctuations. Examples are 
countless. Besides single-cell protein levels [5], they 
range from microbial populations [4] and motor proteins 
with different states strongly bound to the tracks [5] to, 
at a larger scale, insect groups and worms that can switch 
between the directions of motion [6] and, at a still larger 
scale, endemic states of infectious disease, where a part 
of the population is infected while the other part is not 
[7], or stable states in the predator-prey models [5]. 

Even where fluctuations are small on average, occa- 
sionally there occur large fluctuations that lead to dra- 
matic changes in the system. An example is switching 
between dynamically stable states. It plays an impor- 
tant role in biophysical systems including the ones men- 
tioned above. Another example is extinction of a group of 
molecules or a population. A type of extinction that has 
attracted much attention in the literature is spontaneous 
disease extinction in an isolated population, where as a 
result of a fluctuation the number of infected becomes 
equal to zero and the epidemics ceases |M?T]. 

In this paper we will be interested in the control of the 
rates of rare events that result from large fluctuations, 
including interstate switching and extinction. We will 
develop a general approach to control by deterministi- 
cally imposed elementary transitions, in which the num- 
bers of molecules or individuals change. The change is 
small and does not cause switching or extinction on its 



own. We show that, nevertheless, even weak determinis- 
tic control can strongly increase the rates of rare events. 
The approach will then be applied to control of disease 
extinction. Here control is often performed by vaccina- 
tion. Sometimes there is not enough vaccine to eradicate 
the disease. We show how to optimize the control given 
a restriction on the amount of vaccine. The results refer 
also to a more general situation of a limited time-average 
speed at which the transitions are imposed. 

The problem of optimal control of large rare fluctua- 
tions has been discussed in qualitatively different terms 
for two major models of fluctuating systems, continu- 
ous and discrete. Fluctuations in continuous systems 
are often thought of as resulting from an external noise, 
whereas control of such systems is performed by an ap- 
plied regular field. This field can be optimized, for a 
given constraint, so that it will most efficiently change 
the rate of noise-induced rare events |22l |23] . 

Populations or reacting species, on the other hand, 
are intrinsically discrete, with integer numbers of react- 
ing molecules or individuals. The dynamics is often de- 
scribed as resulting from elementary transitions in which 
these numbers change. The transitions happen at ran- 
dom and are characterized by rates. These rates can be 
controlled, but the transitions still remain random. We 
call this elementary transition rates (ETR) control. Opti- 
mal control of this type was discussed in Ref. Because 
of the randomness of the elementary transitions, the con- 
straint on the control field is in some sense probabilistic, 
for example, it can be imposed on the average amount of 
vaccine used in a given time period. 

We will consider here a different type of control, which 
is performed via deterministic processes (DP) where the 
numbers of reacting objects (molecules or individuals) 
change in a well-defined way at well-defined instants of 
time. We call it DP-control. We will show how to in- 
corporate deterministic transitions into the analysis of 
large rare fluctuations. In particular, we will develop 
an eikonal approximation, which applies where the total 
numbers of reacting objects are large. It maps the prob- 
lem of the probabilities of rare events in a reaction system 
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onto Hamiltonian dynamics of an auxiliary system. This 
formulation will be compared with the corresponding for- 
mulation for controlled reaction rates and for controlled 
continuous systems. 

The approach will then be applied to the problem of 
disease extinction with a limited amount of vaccine. We 
show that periodic vaccination implemented in a deter- 
ministic fashion can strongly enhance disease extinction. 
The enhancement can be resonant, for the appropriate 
vaccination period. 



II. MASTER EQUATION IN THE PRESENCE 
OF DETERMINISTIC PROCESSES 



The state of the reaction or population system is de- 
termined by vector X — (Xi,X2, ...) with integer- valued 
components that are equal to the numbers of molecules 
of different sorts or individuals in different population 
groups. The size of the system iV, which gives the typ- 
ical value of |X|, is the large parameter of the theory, 
A'' ^ 1. Our analysis refers to spatially uniform systems, 
which is of interest for many applications in biophysics, 
stirred chemical systems, and also for population dynam- 
ics in moderately large globally connected groups. 
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FIG. 1. Two types of transitions with a change of the num- 
ber of molecules or individuals X — >■ X -I- r. Spontaneous 
transitions happen at random and are characterized by rates 
W(X, r, t). In addition, transitions may be imposed exter- 
nally in a deterministic way. We assume that the duration of 
a transition is small compared to the interval between succes- 
sive transitions, so that the system is described by vector X 
with integer-valued components. 



We will be interested in the probability P(X, t) for the 
system to be in state X. The distribution P changes as 
a result of elementary transitions in which X is changed. 
Such transitions can come from spontaneous chemical or 
biochemical reactions, or from infection and recovery of 
individuals, for example. They can be imposed also in a 
deterministic fashion, see Fig.jl] We will assume that, for 
a given molecule or an individual, the interval between 
successive transitions, whether spontaneous or determin- 
istic, is much longer than the duration of a transition. 
Then elementary transitions can be considered instanta- 
neous, which is consistent with the notion of the compo- 
nents of X being integer- valued. 

In the absence of deterministic transitions the evolu- 
tion of P(X, t) is described by the standard Markov mas- 



ter equation 



where 



5tP(X,t) =T4^P(X,t), 



(1) 



W^P(X, t)=Y. [W^(X - r, r, i)P(X - r, t) 

r 

- W^(X,r,i)P(X,i)]. (2) 

Here, M^(X, r, i) is the rate of an elementary transition 
X — )• X -f r in which the population and/or number of 
molecules changes by r = (ri, r2, . . .). The rates W are 
proportional to the system size, W N . The change of 
X in an elementary transition is independent of N and 
comparatively small, |r| ^ iV, typically |r| ~ 1. We will 
assume that the rates W are either independent of time 
or depend on time periodically. 

We now discuss the effect of a control field that deter- 
ministically imposes elementary transitions. This field is 
a series of pulses applied at instants t\ <ti... < tn < .... 
In each pulse the state vector changes by A, that is, 
X X + A. We will be interested primarily in the case 
where |A| 1. An example is vaccination in popula- 
tion dynamics. One can think that vaccination shots are 
applied at instants t„. In a simplified model where the 
delay between vaccination and acquiring immunity is dis- 
regarded, as a result of a shot the number of susceptible 
individuals decreases by one, whereas the number of vac- 
cinated increases by one. The corresponding components 
of vector A are then -1 and 1, respectively, whereas all 
other components are equal to zero. 

The analysis should be performed differently if a 
change that occurs in a pulse is macroscopic, that is if 
I A| is of order TV. This may be of interest, for example, if 
control is performed by injecting or extracting a macro- 
scopic amount of reacting molecules. The corresponding 
extension is provided at the end of Sec. |III[ 

Disregarding the pulse duration, one can describe the 
effect of the deterministically imposed pulses as a shift 
of the probability distribution. 



P(X + A,t„ + 0) =P(X,t„-0). 



(3) 



The formulation can be immediately generalized also to 
the case where the population change in a pulse A de- 
pends on the instant t„. We will be primarily interested 
in the case where the inter-pulse intervals are larger than 
the typical reciprocal reaction rate. 

In the presence of a pulsed deterministically acting 
field, master equation ([T]) applies in between the pulses. 
Equations (fl} - (pi) fully describe the evolution of the dis- 
tribution P(X, tY However, because of the discontinuity 
of P at instants i„, they are inconvenient for the analysis. 



III. THE EIKONAL APPROXIMATION 

For a large system and | A| ~ 1, the typical inter-pulse 
interval t„ — tn-i is small, cx 1/iV, since the control field 
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is applied to different individuals. If this were not the 
case, the overall effect of the field would be negligibly 
small. On the other hand, of interest is the evolution of 
the distribution P on much longer times, of the order of 
the relaxation time of the system i^. The time charac- 
terizes how the system approaches the stable state in the 
absence of fluctuations and is independent of the system 
size. 

Describing evolution of P on times ^ tj- requires coars- 
ening over time intervals oc and in particular over 
the inter-pulse interval. This can be done by introducing 
the characteristic pulse speed £,{t), 



1 



N{U. 



(4) 



The continuous function ^ {t) is obtained by analytic con- 
tinuation from the discrete set S,{tn)', it is assumed that 
i{t) is smooth, i{t){tn - t„-i) < for |i - t„| < t,. 
Note that we are calling ^(i) speed in contrast to rate to 
emphasize its deterministic nature and to distinguish it 
from the rates used in other equations. 

Function ^(t) determines the number of pulses N^(t)T 
that happen in a time interval r such that t„ — ^n-i ^ 
T ^ tr- We will assume that ^(t) is either constant or 
periodic in time, with period T t„ — We further 

assume, for simplicity, that this period coincides with 
the period of the rates VF(X, r, t) if these rates are time- 
dependent. Then, in the neglect of fluctuations, the sys- 
tem can have periodic steady states; the results can be 
immediately extended to the case where the periods of 
W and ^ are commensurate. 

Introducing ^{t) does not allow one to directly coarsen 
the master equation for the probability distribution 
P(X, t) over time. Indeed, on the tail of the distribution, 
which is of interest for the problem of large rare fluctu- 
ations, function P is steep, |P(X + A,t) - P{X,t)\ ^ 
P(X,<), which means that it significantly changes from 
pulse to pulse and thus -P(X, t + r) — P(X, t) is not given 
by rdtP even for t <^t^. 



A. The Hamiltonian formalism 

A differential equation that describes important fea- 
tures of the evolution of the distribution in coarsened 
time can be obtained by seeking P(X, t) in the eikonal 
form IMS], 

F(X,i) = exp[-A^s(x,t)], xeeX/A/-, (5) 

where x is the scaled state vector. In the limit TV ^ 1 
this vector is quasi-continuous. We assume that function 
s is a smooth function of x and t. This assumption will 
be checked a posteriori. 

To obtain an equation for s between the pulses one can 
write the rate operator in Eq. ([2| as 

WP = [exp(-r(9x) - l]Vt^(X, r, t)P(X, t). 



Functions VF(X, r, i) usually smoothly (polynomially) 
depend on X. In contrast, function P can be exponen- 
tially steep, as seen from Eq. ([s]). Respectively, to the 
leading order in \/N one should only differentiate P in 
the above expression. Then from Eqs. ([T]) and ([s]) one 
obtains 



s(x+^A,t„ + 0) = s(x,t 



0), 

tn < t < tn+1 



(6) 



(7) 



where n = 1,2,..., 

77^ (x, p, t) = V u;(x, r, t) [exp(pr) - 1] , 

^ — 

w{ic,r,t) = j^W{N^,r,t). 

The terms disregarded in the equation for s between the 
pulses give corrections oc 1/N |28j . 

The change of s resulting from a single pulse is small, 
cx 1/N [at the same time, from Eq. ([5]), P can change 
significantly]. This allows us to obtain a single differen- 
tial equation for s on a coarsened time scale. Indeed, 
choosing time r so that tn — in-i r ^ ir, in which 
case many pulses occur during the time r, but the overall 
change of s accumulated between the pulses is small, we 
obtain 

s(x, t + r)- s(x, t) w -THw{yi, 9xs(x, t),t) 

- Tat)A ■ d^s, (8) 

where £,{t) is given by Eq. Q. We can now formally go 
to the limit r — 0, which gives 

9ts(x,t) = -i?„(x,axs(x,t),t)-eWA-axS. (9) 

Equation ([9| has the form of a Hamilton- Jacobi equa- 
tion, with s being the action of an auxiliary dynami- 
cal system with coordinate x, momentum p = d^s, and 
Hamiltonian 



i/(x,p,t) = i/„(x,p,t)+e(i)pA. 



(10) 



Hamiltonian H contains both the well-known part 
that comes from random elementary transitions |26[ 1271 
I29| and a part that comes from the deterministic pro- 
cesses. The latter is described by the last term in H and 
is characterized by two natural parameters: the speed of 
transitions per individual ^(t) and the change in popula- 
tion at a transition A. The structure of this term is very 
different from that of the term H^- it depends linearly 
on the momentum of the auxiliary system p, whereas 
depends on p exponentially, that is, much stronger. Since 
functions w,^ are either independent of time or periodic 
with the same period, Hamiltonian H as a, whole is either 
independent of time or depends on time periodically. 

The reduction of the problem of the distribution tail 
to the Hamilton-Jacobi equation ([£]) makes it possible to 
study the distribution in the presence of deterministic 
control using standard tools of the rare-event theory [26l 
[27| [29] . The assumed smoothness of s as function of 
X immediately follows from Eq. ([9| for smooth initial 
conditions. 
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B. Comparison with the previously studied models 



It is instructive to compare Eqs. ([9| and ( 10 ) with the 
description of large fluctuations in noise-driven continu- 
ous systems in the presence of a control field. There, for 
white Gaussian noise, the probability density can be also 
sought in the eikonal form of Eq. (Is]), with N replaced 
by the reciprocal noise intensity 30j! The Hamiltonian 
has a structure similar to Eqs. ([7| and (10 1, with the dif- 
ference being that exp(pr) in should be expanded to 
second order in p. The form of the control- field induced 
term is the same as in Eq. (10), if the field corresponds 



to an additive force driving the continuous system, with 
-C(i)A playing the role of this force, cf. [HI US]. 

We should also compare the present formulation with 
that for the elementary transition rates control [21] . The 
type of the ETR control most close to the considered here 
DP control is where elementary transitions X — > X + A 
rather than being directly imposed occur at random at 
a controlled rate N^pj.{t) which is independent of X [15] . 
As we will see, in the mean-field approximation where 
fluctuations are disregarded, the evolution of the system 
in the presence of such control is described by the same 
equation as for the DP control. 

To account for the ETR control, one should incor- 
porate the term describing the controlled transitions in 
the operator W in Eq. ([l}, W{X,r,t) W{X,r,t) + 
-^?pr(0'^r A' eikonal approximation, to the lead- 

ing order in 1/A^ one then obtains the Hamilton- Jacobi 
equation for s(x, t) « —N~^\nP{Nx,t), with Hamilto- 
nian 

i7pr(x,p,0 =-ff»(x,p,<)+eprW[cxp(pA)-l]. (11) 



tn+1 — tn '■-^ U. The tail of the distribution can still be 
sought in the eikonal form, Eq. ([5|. However, now the 
equation for s(x, t) is of the form 

dts = ''H^{x.,d:^s,t) 

+ J2Sit-tn) [exp(iV-iA-ax) -1] s. (12) 

n 

It differs from the Hamilton- Jacobi equation, and the 
general analysis of this equation is beyond the scope of 
this paper. However, we are interested here in a com- 
paratively weak deterministic control, where |A|/7V is 
small. Then it is sufhcient to keep the lowest-order term 
in |A|/A^ in Eq. ( [l2| ), and this equation becomes of the 
same form as Eq. , with S,{t) in Eq. (|9| of the form of 
a periodically repeated i5-function. 



IV. MEAN FIELD DYNAMICS AND OPTIMAL 
FLUCTUATIONS 

The action s(x,i), and thus the distribution P to the 
leading order in 1 /N can be found from the Hamiltonian 
equations of motion that follow from Eq. (10), 



x = y rw(x,r,i)eP--+e(t)A, 

P = -V ax«^(x,r,i)(eP'--l). (13) 

Z Jj- 

These equations have an important solution 



ru;(x,r,t) +^(t)A, p = 0. 



(14) 



The major distinction from the DP control, Eq. ( |10[ ), is in 
the different dependence of the control term on momen- 
tum p. The origin of this difference and its consequences 
for the control will be discussed in Sec. |TV| Formally, the 



It describes the mean-field dynamics of the population, 
i.e., the dynamics in the neglect of fluctuations. Equation 
( 14 ) can be obtained also directly from Eqs. ([!]) and ([s]) 
by multiplying them by X, summing up over X while 
disregarding the width of the distribution P(X, t), and 
coarsening over time. In the case of the ETR control 



DP-term in the Hamiltonian Eq. (10), looks like a described by Hamiltonian (11), the equation of motion 



small-pA expansion of the term oc in i/pi- 



C. Simultaneously imposed multiple transitions 

The analysis above should be modified if the control 
pulses are applied rarely but cause macroscopic changes 
in the system, so that a macroscopic portion of the 
molecules or population is changed in a single pulse. For- 
mally, this means that |A| ^ 1 and the ratio |A|/A^ is 
independent of N in the limit of large N. As mentioned 
above, an example is provided by injection or extraction 
of a macroscopic number of molecules of a certain sort 
into the chemical reactor or the biological cell, or a group 
of individuals into the population. If the pulses are short, 
one can describe their effect by Eq. ([3|. We will assume 
that |A|/iV is sufficiently small, so that the system is 
weakly changed by a single pulse. 

For the system to have a steady periodic state, the 
pulses should be applied periodically, with t„ — — 



for X has the same form as Eq. (14), except that £^{t) is 
replaced by ^pi{t). 



We assume that Eq. (14) has a stable solution x^. We 



will consider the simple and fairly common case where x^ 
is stationary, for time-independent rates u'(x, r,t) and 
constant ^(i), whereas for periodically varying in time 
rates and/or periodic ^(i) (with the same period), x^ is 
periodic with the same period as the modulation. The 
time tr gives the relaxation time of x and can be found 
from Eq. ( 14 ) linearized about xa ■ 



We also assume that Eq. ( 14 ) has an unstable saddle- 



type stationary or periodic state X5 on the boundary of 
the basin of attraction to xa- In the problem of extinc- 
tion, at X5 one of the types of molecules or population 
groups becomes extinct. Respectively, one of the com- 
ponents of vector x, which we call xe, becomes equal 
to zero. The state X5 is stable with respect to all com- 
ponents Xi^E- Note that, by construction, Xi > 0, the 
system cannot go beyond the extinction state to nega- 
tive Xe- In addition we assume that, once the group has 
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gone extinct, it does not re-emerge, that is, once the hy- 
perplane xe = has been reached, the system will stay 
there. 

In the problem of switching between coexisting stable 
states, on the other hand, X5 is a conventional saddle 
state. In this case, generally (x^)^ > for all i and the 
system can be on any side of X5 along any component x^. 
From the opposite sides of X5 in the unstable direction 
the system will go to different stable states, in the neglect 
of fluctuations. 

We assume that there is only one state x^ on the 
boundary of the basin of attraction to xa- The system 
can either switch or one group of molecules or popula- 
tion can go extinct. The control field can shift the equi- 
librium position X5. In the extinction problem, the un- 
stable state should be in the hyperplane xe — 0. The 
control field cannot lead to the re-birth of the group that 
has gone extinct. Therefore in the considered model of 
imposed transitions A^; = for the extinction problem. 



A. Logarithmic susceptibility 



Equations ( 13 ) can be used to find the rates of 



fluctuation-induced switching and extinction in biophysi- 
cal and population systems. The problems of calculating 
these rates have much in common, but are not identi- 
cal. The difference between them comes from the fact 
that in switching the system has to cross the boundary 
of the basin of attraction of the initially occupied stable 
state xa and go to the other state, whereas in extinc- 
tion it suffices to reach the hyperplane = 0, and then 
the system will stay in this hyperplane. This difference 
shows in Eqs. (16) and (17 1. 



To logarithmic accuracy, the extinction and switch- 
ing rates We are determined by the probability distri- 
bution P(X,t) for x X5 and t t^, given 
that initially the system was in the vicinity of state x^i 
[UlETj, with X5 being the extinction and saddle states, 
respetively. Taking into account that, in this time range, 
the distribution is maximal at xa , we have [TU [THl [111 HZ] 



M^e ocexp(-Q), Q = iVSext, 

/OO 
dt[p±~ H{x,p,t)] . 
-00 



(15) 



The quantity Sext is given by the action s(x) for x — > X5 
counted off from the value of s at the stable state x^. 
It can be calculated by minimizing the functional s^xt 
with respect to (x(t), p{t)) , which leads to finding the 
Hamiltonian trajectory (13) that starts for t — ^ —00 at 
state XA and approaches state X5 for t — > 00. This opti- 
mal trajectory, (xopt(t), Popt(i))7 gives the most probable 
path followed by the system in spontaneous extinction or 
switching. 

The boundary conditions for the switching and extinc- 
tion problems were discussed in [T51 HH [57]. One can 
extend the analysis to show that, in the presence of the 



DP control that we consider, these boundary conditions 
do not change. Specifically, in the switching problem 

p^O, t^±oo, (16) 

whereas in the extinction problem 

p — 0, t — > —00; Pi^E — > 0, i — 00, (17) 

while the component pe remains nonzero for t ^ 00. In 
the both problems, the Hamiltonian — > for t — > ±00. 

Of utmost interest for us is the effect of weak periodic 
control, 



(18) 



where T is the modulation period. In this case Soxt has 
a simple form Sext — s 



(0) 



.W, where 



.(1) 



m] 



mm 

to 

(0) 



dixit -to)m, 



(19) 



Here, s'^^ is the action in the absence of the control field. 



and (x^pj(i), p^pj(t)) is the most probable (optimal) ex- 
tinction or switching path for ^ = 0; s^^^ gives the field- 
induced correction to the action. We are interested in 
the case where s^xt < oi^ly this case the control 
increases We- 



Equation ( 19 ) is written for the case where the elemen- 



tary transition rates VF(X, r) are independent of time. In 
this case the fluctuation leading to extinction or switch- 
ing can happen at any time with the same probability; 



respectively, the optimal path [xZ{t),p]^J^{t)) can be 

centered (for example, |x;Qp^(t)| can reach maximum) at 
an arbitrary time io- A time-dependent control lifts this 
time degeneracy. It synchronizes the optimal path so as 
to most strongly decrease action Scxt- This synchroniza- 
tion is formally described in Eq. ( 19 ) by minimization 
with respect to Iq. 

Even where the correction to the action is small, 
kcxtl ^ '^ixt: the change of the switching or extinction 
rate, which is given by the factor exp(— iVs^^|.), can be 
very large. It is this large change that makes the con- 
trol exponentially efficient. Function x{t) is called the 
logarithmic susceptibility [3T]. It describes the linear re- 
sponse of the logarithm of the probability P(X) to the 
control field. 



Equation ( 19 ) describes also the effect of simultane- 
ously imposed multiple transitions, where |A| ^ 1, pro- 
vided \p''°Ji{t)A\/N < 1. Function C(t) in Eq. ([l9| 



should be replaced with —N ^ J2n ^{^ ^ ^n) in this case. 

If the elementary transition rates VF(X, r, t) are peri- 
odic in time, the optimal trajectories {yi'^^^{t),p^^^^{t)^ 
are periodically repeated in time, but their phase is fixed 
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by the time dependence of W(X,r,t). If ^{t) has the 
same period, the correction to the action is 

4ii[e(i)] = r dtx{t)m 



with x(t) of the same form as in Eq. (19). In this case 



however, the control field does not synchronize transi- 
tions and its effect depends on its phase with respect to 
the optimal trajectories of extinction or switching. 



B. Comparison with the control of elementary 
transition rates 

The analysis of the weak ETR control was done 
earlier [ISl m] . It is based on the eikonal formulation with 
the Hamiltonian ( [Tl] ) . The control-induced change of the 
extinction or switching rate We is given by Eqs. ( 15 ) and 
(19) in which S^{t) is replaced with ^pr(i) and the loga- 
rithmic susceptibility is replaced with 



XprW = l-cxp[pi°l(t)A]. 



(20) 



By comparing Eqs. ( 19 ) and (l20| one can sec that — Xpr > 
— X, for the same A. This shows that the sensitivity of 
the extinction rate to the ETR control is higher than to 
the DP control. 

The stronger effect of the ETR control for the same 
rate ^pr(i) as the speed ^{t) and for the same change of 
the number of molecules or individuals in a transition A 
is somewhat unexpected. It can be understood qualita- 
tively in the following way. In both cases the control field 
does not cause extinction on its own. It cooperates with 
the fluctuation that leads to extinction. The stronger the 
field the stronger is its effect, but also the effect increases 
with the increasing momentum on the optimal extinction 
trajectory p = (?xS. This momentum gives the steep- 
ness of the probability distribution, as seen from Eq. ^ . 
The steeper the distribution the stronger is the effect of 
changing it by the control field. 

For the ETR control, the control-induced transitions 
have the Poisson distribution. There is a probability that 
within a given time interval r this number will be higher 
than the average number N^p,.{t)T. Where the distribu- 
tion in the absence of control is steep, that is, |p| is large, 
the major effect on We will come from realizations of the 
control with such more frequent transitions. 

In other words, for the ETR control one can think of 
extinction or switching as resulting from two types of 
fluctuations: fluctuations in the absence of the control 
and fluctuations of the control field itself. The extinction 
rate We is determined by the optimal, most probable fluc- 
tuations of the both types. It is this double optimization 
that makes We more sensitive to the ETR control. 



V. OPTIMAL CONTROL BY DETERMINISTIC 
TRANSITIONS: VACCINATION FOR A 
LIMITED VACCINE SUPPLY 

We now consider the problem of optimal control of ex- 
tinction and switching rates by deterministically imposed 
transitions. To be specific, we consider disease extinc- 
tion, with the control performed by vaccination. The 
results are not limited to this model, the only condition 
that we use is that the speed ^{t) > 0. This condition 
follows from Eq. ([3]). We also assume that ^{t) is peri- 
odic. The optimal form of ^(t), i.e., the optimal control 
protocol, depends on the imposed constraints. The con- 
straint that we consider is fairly general, and again, is 
not limited to vaccination only. Essentially, we consider 
a constraint on the number of imposed transitions per 
modulation period. 

As mentioned earlier, a simple model of vaccination is 
where the vaccine is applied to susceptible individuals, 
and in an elementary transition a susceptible becomes 
vaccinated and thus immediately immune to the infec- 
tion. A natural constraint on the speed of vaccination is 
that the total amount of vaccine per period iVS is fixed, 

(21) 



Jo 



where T is the vaccination period, ^{t + T) = £,{t). 

The optimal vaccination speed ^(t) should minimize 
the disease extinction barrier Q subject to constraint 
(21 ). From Eq. (15), optimal ^(t) is given by the solution 



of the variational problem of minimizing the functional 



■Scxt 



where A is the Lagrange multiplier and Scxt[^] is given by 
Eq. ([15|. 

For comparatively weak vaccination, where ^- 



dependence of Soxt is of the form of Eq. (19), minimiza- 



tion with respect to ^(t) can be done in the same way 
as for probabilistic vaccination where the vaccination is 
a Poisson process characterized by rate [24|; such vacci- 
nation is an example of the ETR control. It follows from 



the form of the constraint (21 ) and the condition ^(t) > 



that the optimal vaccination speed as a function of time 
is independent of the logarithmic susceptibility x(t). If 
the system is stationary in the absence of periodic vacci- 
nation, from Eq. ( 19 ) 



> minxT(^) 



dt£(t) = min Yrfi), 

0<t<T 



XT{t)= x{t + nT). 



(23) 



The minimum of s^^|. is reached, and the inequality (|23|) 



becomes an equality, for ^(t) of the form of periodically 
repeated 5-pulses, 



nT). 



(24) 
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The instant t^in determines where vaccination pulses are 
apphed. It is arbitrary for a stationary system. In fact, 
the disease extinction events are adjusted to the vacci- 
nation pulses. This adjustment is described by the min- 
imization over t in Eq. ( |23[ ); the periodic function XT(i) 
is minimal at imin- 

The optimal shape of ^(t) is easy to understand: the 
vaccine is most efficient if it acts where |xT(i)| is maxi- 
mal, and also where xt(^) < to assure that s^^t < 0. 
We note that strong periodic vaccination has been inves- 
tigated in the framework of deterministic epidemic mod- 
els, where all fluctuations are disregarded, and it was 
found that pulsed vaccination is advantageous compared 
to vaccination at a constant rate |32) . 

For periodically modulated systems, where W^(X, r, t) 
are periodic functions of time, optimal ^{t) still has the 
form of pulses, Eq. ( 24 1 , but now spontaneous extinc- 
tion events are synchronized without vaccination, and 
this is the vaccination that must adjust in order to in- 
crease the rate of disease extinction. Equation (23) for 



Sgxt applies only if tniin is chosen so that xt(0 is mini- 
mal at i niin 5 

a wrong choice of tmin will be less efficient 
and can even slow down spontaneous disease extinction 
instead of speeding it up. 

It is instructive to compare the results with the prob- 
abilistic ETR-type vaccination, where vaccination is ap- 
plied at random with rate S,pi{t). The differences are in 
the form of the logarithmic susceptibility, cf. Eqs. (19) 
and (20), and in the meaning of the parameter S. In 
the deterministic scenario, iVTS gives the actual num- 
ber of individuals vaccinated in time T, which is fixed 
and is determined, for example, by the periodically sup- 
plied vaccine. In contrast, in the probabilistic scenario 
NTEpr = N dt^pr{t) gives the average number of vac- 
cinated individuals per time T. The actual number is 
random. For large NT'E.p,-, the distribution of the number 
of vaccinated is close to Gaussian, with variance TVTSpr. 



A. Resonant optimal vaccination 



Expressions (15), (19), and (24) allow one to investi- 
gate the effect of optimal deterministic vaccination on the 
disease extinction rate, to study how this effect depends 
on the interrelation between the parameters of the sys- 
tem and vaccination, and to find the change of the disease 
extinction rate for specific models. The analysis is simi- 
lar to that for probabilistic vaccination and we will 
not reproduce it here. Generally, the effect increases with 
the increasing period T, as seen from Eq. (23). However, 



there may be important features that require special at- 
tention. 

As an illustration we show the results for determinis- 
tic vaccination in one of the broadly used epidemic mod- 
els, the susceptible- vaccinated-infected-recovered (SVIR) 
model [121 1321 ■ Here, in the absence of vaccination, 
susceptible individuals, with population Xi = S, are 



brought in at rate fiN (birth) , each population decreases 
at rate /iX^, i = 1,...,4 (dying), the infection rate is 
/3X1X2/N, where X2 = / is the number of infected in- 
dividuals, and infected can recover at rate 7X2. The 
vaccination that we discuss corresponds to the decrease 
of the number of susceptible individuals at speed N£^{t). 

We assume that both the recovered (R) and vacci- 
nated (V) individuals keep the immunity, they do not 
become susceptible again. These groups of population 
are "sinks", there is no influx to other groups and the 
transition rates are independent of R and V. The dy- 
namics is then determined by two variables, Xi — S and 
X2 — I. In the mean-field description, for /3 > 7 -f /i 
the model possesses a single endemic state x^; for fi < 
4 (/3 — 7 — /i) (7 -t- /i)^/3~^ this state is a focus on the plane 

{xi,X2). 
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FIG. 2. (a) The Fourier transform of the logarithmic sus- 
ceptibility with respect to deterministic vaccination for the 
SVIR model, Eq. (25 1. The parameters are 13/ jj. — 80 and 



■y/^j, — 50. The susceptibility spectrum displays a sharp peak 
at the characteristic frequency ljo of decaying oscillations near 
the mean-field endemic state, (b) The change of the scaled ex- 
tinction barrier s' = /xSg^t/— with vaccination period T. The 
solid line shows s' where there is no limit on vaccine accumu- 
lation. The dashed lines refer to the case of a limited amount 
of accumulated vaccine M, where the actual mean speed of 
vaccination is Ha = min(H, M/T). The scaled accumulation 
limit is M' = fiM/E,. The locations of resonances of s' are 
independent of M. 

In the mean-field approximation, the populations of 
susceptibles and infected exhibit decaying oscillations in 
time as the system approaches the endemic state. In 
the same parameter range, the populations oscillate also 
on the optimal disease extinction path [H]. Because of 
the oscillations, the Fourier-transform of the logarithmic 
susceptibility 



dtx{t) exp(za;i) 



(25) 



displays a resonant peak shown in Fig. [2] (a) [x(i^) is 

obtained using the optimal trajectory ^XQpt(t), pQpt(t)^ 

found in Ref.|24]. The peak of x{^) is centered at the fre- 
quency ujq of decaying (in the mean-field approximation) 
oscillations near :x.a- 
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Because of the spectral peak in x(ijj), the vaccination- 



induced term in the extinction exponent Q^^' = -/Vs^^j). 
depends on the vaccination period T nonmonotonically. 
This behavior is illustrated in Fig. [2] (b). Where T is 
close to a multiple of the oscillation period 2tt /ujq, "sixt 
displays a peak. Respectively, the disease extinction rate 
is exponentially enhanced in this case. 

Shown in Fig. [2] (b) are also the results for s^^2t in 
case where the total amount of accumulated vaccine is 
limited, which means ST < M, where M characterizes 
the accumulation limit. Such constraint is typical for live 
vaccine, as it may be dangerous to store too much vaccine 
in this case. For a given period T, this is effectively 
a constraint on the average vaccination speed 2, as it 
makes no sense to increase it beyond M/T. If we set 
S = M/T for large T, then it is seen from Eq. ([23]) that 



the effect of vaccination saturates with increasing period. 



.(1) 



Af minx(t), 



However, the maximum of |Soxtl reached not for T — > 
00, as seen from the dashed lines in Fig. |2](b), but for T 
close to an appropriate multiple of I-k jujQ. 

If the maximum of Isp^t I li^^ where the solid and dashed 
lines are separated, as for M' = 3, the saturation limit 
has already been reached for the corresponding T, and 
the actual vaccination speed Sa = min(S, M/T) is equal 
to M/T. On the other hand, for M' = 1 the maximum 

of |sp^j| lies practically on the solid line, and thus it is 
reached for « S. In practice, for a given maximum 
accumulation level M and a given range of available av- 
erage speeds S and vaccination periods T, one should 
adjust S and T so as to take advantage of the resonance 
while keeping Sq equal to S. 



T 



00. 



VI. CONCLUSIONS 

We have developed a theory of control of large rare 
fluctuations in reaction and population systems by de- 
terministically imposed transitions. These transitions 
occur at well-defined instants of time and cause well- 
defined changes in the numbers of individuals or reacting 
molecules. The control is comparatively weak, so that 
the mean-field dynamics of the system is perturbed very 
weakly. Nevertheless the probabilities of rare events may 
be changed significantly. 

The theory applies to mesoscopic systems, where the 
characteristic number of molecules or individuals N is 
large, but not exceedingly large, so that the rates of rare 
events, which depend on N exponentially, are not ex- 
ceedingly small. The analysis takes into account the dis- 



creteness of the numbers of individuals or molecules and 
the fact that the tail of the probability distribution is not 
smooth, the distribution changes by a factor where the 
number of individuals/molecules changes by one. Our 
approach is based on the eikonal approximation. It al- 
lowed us to reduce the problem of the distribution tail to 
a Hamilton- Jacobi equation for the action of an auxiliary 
conservative dynamical system and express the Hamilto- 
nian in terms of the characteristic speed of the determin- 
istic transitions. 

Of primary interest for this paper was control of such 
rare events as switching between coexisting stable states 
or population extinction. We found that even a com- 
paratively weak control field can lead to an exponential 
increase of the rates of these events. The exponent is 
proportional to N and linearly depends on the speed of 
the imposed transitions. It is determined by the motion 
of the system in the most probable fiuctuation leading to 
the event in the absence of the control. 

The considered control should be contrasted with the 
ETR control, where all elementary transitions happen at 
random, but their rates can be controlled. Unexpectedly, 
we found that the ETR control can be more efficient than 
the control by deterministically imposed transitions, pro- 
vided the average transition speed and the composition 
change in a transition are the same. This is a consequence 
of the double optimization in the ETR control, where the 
rate of the rare event is determined by both the optimal 
fiuctuation leading to this event in the absence of control 
and the optimal realization of the control. 

We have considered the problem of optimal control of 
rare events where a number of deterministically imposed 
transitions per period is limited. The specific example 
is disease extinction by vaccination in the situation of 
a limited average speed of vaccine supply. The optimal 
shape of the vaccine pulses is a train of (5-pulses. It is 
independent of the model of the system. The results are 
illustrated using the well-known susceptible-vaccinated- 
infected-recovered model, and the possibility of resonant 
exponential enhancement of the effect of vaccination is 
demonstrated for the deterministic vaccination. 

We note in conclusion that the theory of rare events 
in mesoscopic population and reaction systems developed 
in this paper fills the gap between the previously studied 
deterministic control of rare events in dynamical systems, 
which are continuous, and the ETR control of reaction 
systems, which are inherently discrete. 
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